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We study the excess free energy due to phase coexistence of fluids by Monte Carlo simulations 
using successive umbrella sampling in finite LxLxL boxes with periodic boundary conditions. Both 
the vapor-liquid phase coexistence of a simple Lennard- Jones fluid and the coexistence between A- 
rich and B-rich phases of a symmetric binary (AB) Lennard- Jones mixture are studied, varying the 
density p in the simple fluid or the relative concentration xa of A in the binary mixture, respectively. 
The character of phase coexistence changes from a spherical droplet (or bubble) of the minority phase 
(near the coexistence curve) to a cylindrical droplet (or bubble) and finally (in the center of the 
miscibility gap) to a slab-like configuration of two parallel flat interfaces. Extending the analysis of 
M. Schrader, P. Virnau, and K. Binder [Phys. Rev. E79, 061104 (2009)], we extract the surface 
free energy 'y(R) of both spherical and cylindrical droplets and bubbles in the vapor-liquid case, 
and present evidence that for R — > oo the leading order (Tolman) correction for droplets has sign 
opposite to the case of bubbles, consistent with the Tolman length being independent on the sign 
of curvature. For the symmetric binary mixture the expected non-existence of the Tolman length 
is confirmed. In all cases and for a range of radii R relevant for nucleation theory, 'y(R) deviates 
strongly from 7(00) which can be accounted for by a term of order j(oo)/y(R) — 1 oc RT 2 . Our 
results for the simple Lennard- Jones fluid are also compared to results from density functional theory 
and we find qualitative agreement in the behavior of j(R) as well as in the sign and magnitude of 
the Tolman length. 

I. INTRODUCTION 

Curved interfaces between coexisting vapor and liquid phases (or between coexisting A-rich and B-rich phases of 
binary (A,B) liquid mixtures having a miscibility gap) are ubiquitous in nature. Nanoscopic spherical droplets or 
bubbles need to be considered in the context of nucleation phenomena [H-Q • In nanoscopic slit pores or cylindrical 
pores, wetting or drying phenomena [H-Q typically can cause a curvature of interfaces that extend across the pore 
[Tol - fT^ I . Note that phase coexistence in porous media is relevant for widespread applications [IB4I3 ; ranging 
from oil recovery out of porous rocks to devices in microfluidics. 

Understanding the properties of such curved interfaces is a longstanding and difficult problem. Note that on 
the atomistic scale interfaces between coexisting phases are rather diffuse and the problem of understanding their 
structure is not fully solved. From the point of view of thermodynamics, the central quantity of interest is the 
intcrfacial tension 7 and its dependence on the radius of curvature R of the droplet (or bubble). Although this 
problem was already mentioned by Gibbs [HI and discussed in classical papers by Tolman [l9[ , the subject is still 
controversial. Tolman introduced [l9[ a length S(R), referred to as "Tolman's length" , to describe the curvature 
dependence of j{R), but the understanding of this length is incomplete till now [20l - l42l ]. 

When one considers a droplet coexisting with surrounding vapor, different definitions of the droplet radius are 
conceivable. One of them is the "equimolar radius" R, which assumes that the volumes of the two phases Vj , Vn 
and their particle numbers Ni, Nn arc additive, 

V = Vr + V n , N = Ni + N11, (1) 

such that there is neither an excess volume nor an excess particle number associated with the interface. This 
definition implies that the "intcrfacial adsorption" is identically zero and in equilibrium both phases / (liquid) and 
II (vapor) have the same chemical potential /!/ = p,n = /i, since they can exchange particles. For a spherical 
droplet (or bubble) we have Vj = 47ri? 3 /3, of course, and in the grand-canonical ensemble with fj, and temperature 
T as control variables, the densities of the coexisting phases pi(/j,,T) = N1/V1, pn(p,,T) = Nu/Vu are those 
of bulk liquid in equilibrium and surrounding metastable vapor. Here we disregard the conceptual problem that 
metastable states in the framework of statistical mechanics arc not completely well-defined 0, |43[ . 

Another important concept to introduce a division between the two phases is the "surface of tension", i.e., the 
surface where the surface tension acts fl8l - [2l"l ]. Consider, for fixed R as defined above, the surface tension j(R p ) 
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[18H21L l43l|. i.e., the excess contribution to the thermodynamic potential, as a function of the radius R p where 
we put the dividing surface. One expects that j(Rp) will exhibit a minimum for a choice of R p somewhere in 
the region of the atomistically diffuse interface, but in general there is no principle that requires that R p and R 
exactly coincide. Considering the pressure difference Ap between the coexisting phases, one can derive a generalized 
Laplace equation [20| 

Ap = 2 1 {R p )/R p + d 1 {R p )/dR p . (2) 

Thus at the "surface of tension" where d^(R p )/dR p = the pressure difference takes its standard macroscopic 
form, Ap = 2j/Rp. 

Now the proper definition of the Tolman length can be written as 

5= lira 5(R)= lim (R-Rp). (3) 

R— >-oo R— foo 

Recall that in this treatment R p and hence S(R) in general can depend on R. Tolman [l9| also suggested the 
approximation that 6 (for a liquid droplet) is a positive constant. If this is assumed, one can show that the 
curvature-dependent surface tension becomes 

7( J R)= 7 (oo)/(l + 25/i?). (4) 

On the other hand, there is evidence from density functional calculations [2!| [3(| [H, H3 for liquid droplets 
surrounded by supersaturated vapor that actually S(R) varies strongly with R, even changing its si gn from a 
positive value at small R to a small negative value at large R. Apart from very recent indications fill |42| , most 
simulations, some of which are still inconclusive (see the discussion in |I(J), did not support this result. (Note that 
the positive values for 6 in some previous simulations have not been obtained by usin g: Eg. (EH but by a "virial 
route" through l/i?-corrections to the pair correlation function of a planar interface [23, |4(|-) There are also 
compelling arguments (22l. l26j that for systems obeying a strict symmetry between the two coexisting phases, e.g., 
the Ising lattice gas model of a fluid that exhibits particle/hole symmetry, the Tolman length S must be identically 
zero, since interchanging the identity of the coexisting phases turns a droplet into a bubble which simply means 
a change of sign of the radius of curvature of the interface separating them. This implies that j(R) for such 
symmetric systems can only be a function of R 2 , invalidating Eq. ((4]) even for arbitrarily large R. 

There appears to be consensus about the uniqueness of the 1/7? form of the leading correction to the surface 
tension of bubbles and droplets (given by the Tolman length) when analyzed in different frameworks. This is 
not so for the next-to-leading correction and even for the leading correction to the surface tension of a cylindrical 
interface [26|,[44j]. Nevertheless, a phenomenological expansion of the surface tension of an arbitrarily curved surface 
in powers of its curvatures can be devised using the Helfrich form. For spherical (subscript s) and cylindrical 
interfaces (subscript c) the following expressions are obtained [3l| 

ls (R) = 7(00) - 27(00)- + (2k + k)-L , (5) 
6 k 1 

lc(R) = 7(00) _ ^(°°)^ + 2^2 • ( 6 ) 

Here, h is the bending rigidity constant and h is the rigidity constants associated with Gaussian curvature. Such a 
form neglects possible nonanalytic terms in R, and its usefulness should be judged by comparison to actual results. 

In the present work, we make an attempt to study the problem of the Tolman length and of higher order 
corrections, both by analyzing computer simulations of simple models for fluids and fluid mixtures, and by density 
functional theory. The distinctive features of our work are that we apply a recent "thermodynamic" method [4JJ 
based on exploiting Eq. (JTJ) for various finite volumes V = L 3 with periodic boundary conditions throughout and 
calculate "f(R) directly via application of successive umbrella sampling methods [45[ to obtain accurate estimates 
for the appropriate thermodynamic potential of our model systems. In subsequent sections we use the symbols 
jab for surface tension along the A — B interface in a binary mixture and j v i for vapor-liquid interfacial tension 
in a single component system. 

The paper is organized in the following sequence. To clarify the curvature dependence in symmetrical situations, 
in Sec. II we present results from the study of a symmetric binary (A,B) Lennard Jones mixture whose equilibrium 
properties have been extensively studied in earlier works [46|, [47]]. In Sec. Ill, we turn attention to the vapor-to- 
liquid transition of the Lennard- Jones fluid, considering both droplets and bubbles on an equal footing. Previous 
works, with other methods, have considered droplets almost exclusively. However, our simulation setup also allows 
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us to consider cylindrical droplets or bubbles, stabilized by the periodic boundary conditions. Such cylindrical 
interfaces are not only of interest to provide further constraints on the possible value of S as defined in Eq. ([3]), 
but are also important when one considers phase coexistence in slit pores. Sec. IV discusses the problems from 
the point of view of density functional theory, while Sec. V summarizes our conclusions. 

II. THE CURVATURE-DEPENDENT SURFACE TENSION IN A SYMMETRIC BINARY 

LENNARD-JONES MIXTURE 

We study a binary fluid of N point particles (labeled by index i, at positions ft in the cubical box of finite volume 
V = L 3 subject to periodic boundary conditions) interacting with pairwise potentials ufaj) with = |r, — rj\. 
Starting from a full Lcnnard- Jones potential 

4>Lj(nj) = 4£ Q ^[(cr Q(3 /ry) 12 ~ (^Mj) 6 ], a,P € A,B, (7) 
we construct a truncated potential as follows [47| : 

u{nj) = 4>Lj{nj) - (t>Lj{r c ) - (nj - r c ) ^ LJ | ry=rc ,for nj < r c , (8) 

CLTij 

while u(rij > r c ) = 0. This form ensures that both the potential and the force are continuous at r = r c (48|. The 
potential parameters are chosen as 

paa = °bb = o A B = c, r c = 2.5a; (9) 
saa = sbb = 2e AB = e. (10) 

Choosing a reduced density p* = 1, where 

p* = pa 3 = Na 3 /V, (11) 

we work with a dense fluid and at the temperatures of interest neither the vapor-liquid transition nor crystallization 
is a problem. The unit of temperature T is chosen such that e/ks = 1, also throughout the paper we set the unit of 
len gth a to unity. The critical temperature T c of phase separation into an A-rich phase and a B-rich phase occurs 
at [47| T = 1.4230 ± 0.0005. Using Monte Carlo methods in the semi-grandcanonical ensemble and applying a 
finite-size scaling analysis [H, H^, the phase diagram in the plane of variables (T, xa = Na/N) has been obtained 
rather accurately in the critical region [47J and we have extended it here to lower temperatures as depicted in FigfT] 
The first task then is to obtain reliable estimates for -fAB = Iab(oo), the interfacial tension between flat, 
infinitely extended coexisting A-rich and B-rich phases. Here we follow the standard method [5ll - [59l | to first sample 
the effective free energy V/l(xa,T) — —kBTh^P/^^VTixA) / Pa^nvt{x c a°^)\, presented in FigEl by successive 
umbrella sampling in the semi-grandcanonical ensemble where the chemical potential difference Ap between A and 
B particles is independently chosen. Note that the flat region (hump) of Jl (xa , T) in Fig. [2] can be interpreted as 
being caused by the excess free energy of two interfaces, of area I? each, oriented parallel to two axes of the cubic 
box, separating the A-rich domain from the B-rich domain, and hence 

7 abM = lim jab(L), lAB {L) = Lf L (x A » 0.5)/2. (12) 

Note that in finite boxes the capillary wave spectrum [B-d, [2(| of the interfaces is truncated, and there may be 
additional corrections due to the translational entropy of the interfaces, residual interactions between them, etc. 
Therefore, an extrapolation of "/ab(L) to L — > oo is indeed important to reach a meaningful accuracy. Such exercise 
is shown in FigJ3] where "/ab{oo) = 0.722 is obtained from a linear extrapolation as a function of inverse linear 
dimension L of systems. 

As described in the work of Schrader et al. [4l| and Winter et al. [5^, [6(|, the curvature dependant surface 
tension "] A b{R) is extracted from /l(xa, T) using other parts of the same curve in Fig. [21 namely parts that reflect 
the coexistence of an A-rich droplet with B-rich background. Such states are found in the ascending, size-dependent 
part of fL^XA-,), seen in Fig. 01 and correspond to spherical droplets for smaller xa and cylindrical droplets for 
larger x A , representative snapshots of which are shown in Fig. O In order to be able to analyze /l(xa,T) and 
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distinguish what fraction of this effective free energy is due to the surface free energy of the droplet and what is 
due to the background, a B-rich phase supersaturated in A-particles analogous to a vapor phase in the gas-liquid 
context, we also need to study the effective chemical potential difference A(j,l(x a ): 

-±-A» L (x A ) = [df L (x A ,T)/dx A ] T , (13) 

Kb-1 

which we present in Fig. [6] One can clearly identify the various transitions: the peak in the A[Il{x a ) vs. x A curve 
for small x A is the transition from the homogeneous state of the box to droplet + "vapor" coexistence (the so-called 
droplet evaporation/condensation transition fill I59U641 ]): the next step (near x A ~ 0.2) signifies the transition of 
the droplet from spherical to cylindrical shape; and finally, near x A ~ 0.35, the transition from the cylinder to the 
slab configuration occurs for which we have Ahl(x a ) = 0, of course. In order to extract information on droplet 
surface free energies only those parts of the Jl{x a , T) vs. x a and A^t, (xa : T) vs. x a curves can be used that are 
not at all affected by fluctuations associated with these transitions [4l|, |51h62| . 

Fig. \7\ now recalls the method by which the surface free energy F s (R) of the droplet and its radius R are extracted 
fill l59U6Cl| . An A-rich droplet coexists with a B-rich background having the same chemical potential as another 
state of pure B-rich phase. The chemical potential of the latter must be equal to those of the A-rich droplet. Use 
of Eq.([T]) implies that one takes the bulk properties of the droplet, as a definition, equal to those of this B-rich 
phase, including its bulk free energy density fo as indicated in the lower part of the figure. The difference Af in 
Fig J7] simply is due to the surface free energy of the droplet 

VAf = A-wR 2 lAB {R). (14) 

On the other hand, from Eq. ([1}, or its analog for a binary mixture, we can readily obtain R from the concentration 
difference Ax between the respective states as 

Ax = (l-2xX° x )(47ri? 3 /3£ 3 ), (15) 

where x™ ex is the concentration of the pure B-rich phase on the coexistence curve. Of course, for more precise 
estimation, 1 — 2a;™ cx in Eq. (|T51) should be replaced by the difference between the compositions of the pure A-rich 
and B-rich phases at the considered value of Afj,, but in practice Eq. (|1 5|) is an excellent approximation. Estimating 
hence Af and Ax from the simulation, Eqs. (fl~4|) . (fT5|) yield the corresponding values 7ab(-R) and R. One clearly 
sees that both AfiL(x A ,T) and fL(x A ,T) depend on the linear dimension L of the system, which also appears 
explicitly in Eq. (|15[) . However, the interpretation in terms of ~/ A b{R) makes only sense if this dependence on L 
completely drops out - at least to a very good approximation. 

Fig. IS presents the resulting plot of Fs{R) = 4:TtR 2 j A b(R) vs. R and compares it with the capillarity approxi- 
mation, Fg(R) = 4:TtR 2 j A b(oo). One sees that the latter is accurate when the resulting surface free energies are 
rather large, up to 600 fc^T! According to the classical nucleation theory (CNT) (H-Q, the resulting barrier F* to 
be overcome in a homogeneous nucleation event would be F* = Fs(R*)/3, R* being the radius of a critical nucleus. 
This implies that for 4 < R* < 8 barriers F* /ksT of about 30 to 200 result. On the other hand, one can already 
see from Fig. however, that for very small nuclei the relative deviation from CNT increases. 

We now turn to the question, whether this deviation can be accounted for in terms of Tolman's hypothesis, 
Eg. ((4]). To test for the latter, one notes that a plot of ^(oo)/^f(R) should be equal to 1 + 25/ R, i.e. a straight line 
when plotted against 1/R. However, using the data from Fig. [8] no straight line vs. 1/R can be observed, while 
the data are compatible with a straight line when plotted against 1/R 2 , as seen in Fig. [S] This result agrees with 
the general expectation, discussed already in the introduction, that for systems that exhibit a precise symmetry 
between the coexisting phases the Tolman length as defined in Eq. §5§ is exactly zero [22|, [26| . Our results shown in 
Fig. IH1 as well as a related study for the Ising model implies that instead of Eq. (Q| the curvature-dependence 
of the surface tension 7ab(-R) for symmetric mixtures can be described by 

1ab{R) = 7abM/[1 + 2(£ S /R)% (16) 

where the characteristic length I s at low temperatures is close to the LJ diameter a. It remains a challenge for the 
future to clarify the temperature dependence of this length as T — > T c , however. We speculate that in this limit 
l s simply becomes proportional to the correlation length £ of order parameter fluctuations. This is supported by 
explicit mean-field results for the rigidities k and k [65[ to which £ s is related through C 2 = — (k + k/2)/^j A B{oo) 
(see Eq. ©) . 

In Fig. [5] deliberately only the range 1/R 2 < 0.12 is shown. It should be noted that our method to construct 
7ab(-R) is well-defined in the limit R — > oo, but as R gets small, systematic errors arise: one reason is that 
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there is a nonzero, albeit very small, probability that a state containing a droplet (or bubble), such as shown in 
Fig. [5^, undergoes a fluctuation to cylindrical shape (Fig. [5p) or to a homogeneous phase (Fig. upper part, 
left most cartoon). These rare fluctuations are completely negligible for large L used in Fig. [9] but should become 
important if linear dimensions of the order of a few a were used, and such linear dimensions would be needed 
if we were to continue the study of the behavior towards larger l/R 2 . In addition, artefacts due to the periodic 
boundary condition arc to be expected, when L is only of the order of a few a: fluctuations on the right side of 
the droplet would interact with fluctuations on the left side of its periodic image. Note that droplets with small 
R require the use of small L, in order to ensure their stability. In this respect, our method that does not restrict 
statistical fluctuations beyond the use of periodic boundary conditions on the scale L, differs from the density 
functional theory (DFT) where one can consider the equilibrium of an arbitrarily small droplet (at the top F* of 
the nucleation barrier) in constrained equilibrium with surrounding mctastablc phase extending infinitely far away 
from the droplet. Due to the mean field character of DFT, this constrained equilibrium does not decay while a 
corresponding computer simulation clearly would result in an unstable decaying situation. 



III. DROPLETS VS. BUBBLES IN THE SINGLE-COMPONENT LENNARD-JONES FLUID 



In this section, we focus on the vapor-liquid transition of simple fluids using the LJ potential, following the work 
of Schrader et al. [4l|. But here, in addition to spherical droplets, we also study spherical vapor bubbles as well as 
droplets and bubbles of cylindrical shape. While simulations of droplets surrounded by vapor are truly abundant 
(and such simulations have been attempted since decades 0), other geometries have found little attention, so far, 
despite their phy sical significance. 



As in Ref. [41| we use a simple truncated LJ potential, 



tt(r i3 -) = AeKa/n^-ia/njf+C], r < r c = 2.2 1 / 6 ( r, 

u(nj) = 0,r>r c , (7= 127/16384. (17) 

Thus, the potential is cut at twice the distance of the minimum, and the constant C is chosen such that u(rij) 
is continuous for r%j = r c . Two temperatures T = 0.68T C and T = 0.78T C (T c = 0.999) were studied with the 
linear dimension of the simulation box varying from L = 11.3 to 22.5. The accuracy of the estimation of the excess 
free energy hump /l(p/T), p (= N/V) now being the particle density, for large L is enhanced by applying the 
Wang-Landau method [66j in addition to successive umbrella sampling. Simulation method and the data analysis 
[4~fl ] are rather similar to the description provided in the previous section. 

As an example, from the "raw data" obtained from these simulations, Figs . [TUl and 1111 show the free energy hump 
Af(p,T)/ksT at the chemical potential p = p cocx yielding phase coexistence between uniform saturated vapor 
and liquid, as well as the derivative pl(p) /ksT = (<9(/l(p, T)/kBT)/dp)T, f° r T = 0.78T C . Similar to the binary 
LJ mixture, the flat parts of Jl(p) for p around 0.35, where Ap p(p) = 0, can be used to find the vapor-liquid 
interfacial free energies 7„i(oo) by an extrapolation to L — > oo [5ll - l57| . in full analogy with Eq. (TT2"j) . Fig. [T2] 
presents the counterpart of Fig. [3] in the present case, for T = 0.78T C . 

The data shown in Figs. 1101 and [TT1 are analyzed in an analogous manner as explained in the context of Fig. [7] 
The extension of this method to the case of bubbles is completely straightforward: one simply uses the part of 
the pl(p) and Jl(p) curves near p = 0.6 rather than near p = 0.1 in Figs. [TUl EU While traditionally droplets 
have been identified as clusters of connected particles (e.g. [32[), such methods are not at all straightforward to 
generalize in order to identify bubbles: the present thermodynamic method clearly has an advantage here. Also the 
extension to cylindrical surfaces is straightforward - we simply have to replace Eq. (|14[) by an expression involving 
the surface area 2ttRL of a cylinder surface, 

VAf^ = 2nRL 7 %\R), (18) 

and Eq. (|15p becomes modified by an expression containing the volume of the cylinder, irR 2 L, rather than that of 
the sphere, AttR 3 /3, so that 

Ax = (pi - Pv )ttR 2 /L 2 . (19) 

Fig. E2 presents the counterpart of Fig. [5]for the single-component fluid, plotting Fs/ksT vs. the sphere radius, 
both for droplets and bubbles, and compares them to the capillarity approximation (CNT) where, of course, there 
is no difference between droplets and bubbles. Indeed, we see that CNT overestimates the correct surface free 
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energies somewhat, and Fs(R) for bubbles falls clearly below the result for droplets (such a difference cannot occur 
in symmetric models). Fig. [14] presents the results for cylindrical geometries. 

Figs. H"5l and IT6l show our attempts to extract a Tolman length from the results of r y v i(R) for both droplets and 
bubbles with spherical as well as cylindrical shapes. Motivated by the result for symmetric systems, Eq. (fl"6|). 
where a quadratic term in 1/R is present while the linear term being absent for symmetry reasons, we now assume 
the presence of both linear and quadratic terms and attempt to fit data to the forms: 

Jvl(oo)/y v i{R) = l + 25/R + 2(£ s /R) 2 , spherical droplet, 

Jvl(<x>)hvl(R) = 1 -25/R + 2(£ s /R) 2 , spherical bubble, 

yvl(°°)/lvl(R) = l + 5/R + 2(£ c /R) 2 , cylindrical droplet, 

j v i{oo)/jui{R) = 1 -6/R + 2(£ c /R) 2 , cylindrical bubble. (20) 

Eqs. (|20[) assume that in the leading order, droplets and bubbles, where one goes from a convex to a concave 
interface, just differ by a change of sign. The lengths £ s and £ c are related to the rigidities k and k introduced in 
Eqs. © and © by 

,2 _ 25 2 _ 2k + k 

S 2 k 

11 = y - ' (22) 

Upon comparing Eqs. (|20p to the actual results obtained from unconstrained fits to the simulation data, quoted 
in Figs. 1151 and [TBI we observe a clear indication of linear correction, having a positive sign for bubble while being 
negative for droplet. This tentatively implies a negative Tolman length of order 5 « —0.07 ± 0.04 for T = 0.68T C 
and 5 » —0.11 ± 0.06 for T = 0.78T C . Of course, the large error bars which we simply extract from the scatter of the 
four individual estimates (droplets and bubbles of spherical and cylindrical shapes) assuming that the symmetries 
postulated in Eqs. (|20|) holds, are somewhat disappointing. Of course, the possibility of systematic effects due to 
higher order terms oc i?~ 3 ,i?~ 4 etc. cannot be ruled out. However, the success of Eq. (|16p in the symmetric case 
for a similar range of 1/R (Fig. [9]) strengthens our expectation that for the present case these higher order terms 
are still negligible. The coefficient of the term 1/R 2 is of order unity, i.e., the lengths £ s and £ c are of the order of a. 
as found in the symmetric case. From our fits we find negative bending rigidities k with a magnitude of about half 
a ksT whereas the rigidity k is consistent with zero. We also note that the lengths 5 and £ s , £ c increase in absolute 
magnitude, with increasing temperature. Thus, at least there is no qualitative contradiction with the prediction 
that actually 5 should diverge as T — > T c (39|. However, the smallness of S at T = 0.78T C precludes any hope that 
this possible divergence might be probed by simulations. 

We also note that ten Wolde and Frenkel [12] in the analysis of their simulation results for droplets (applying a 
rather different method than in the present paper) sug gested that j v i(oo)/y v i(R) — 1 oc 1/R 2 . This was motivated 
by theoretical results of McGraw and Laaksonen [68|, assuming hence that 6 = for a Lennard- Jones fluid. In 
view of our result, that 5 clearly is an order of magnitude smaller than the length £ s controlling the magnitude 
of the R~ 2 term, it is understandable why ten Wolde and Frenkel missed the existence of the Tolman correction. 
The strong point of the present work is the proof of a clear difference between 7„;(-R) for droplets and bubbles, 
which is inconsistent with [f38j]. Figs. I13I14I15I16I give very clear evidence for the existence of this difference, which 
according to Eq. ([2H|) should be 

A7 = [7«i(oo)/7„z(i?)]bubbics - [lvi{oo)/j v i(R)] dmp i cts — AS/R + 0(R~ 3 ), spheres 
A 7 = [7w(°o)AM(-R)]bubbics - [-fvi{oo)/-fvi(R)}di- op icts = 2S/R + 0(R~ 3 ), cylinders (23) 
Our numerical data would yield, for T = 0.78T C , 

A7 = 0.3811(l/i?) + 0.266(l/i?) 2 , spheres, 

A 7 = 0.2386(l/i?) - 0.172(l/i?) 2 , cylinders. (24) 

The small value of the coefficient of 1/R 2 , compared to the individual contributions coming from droplets and 
bubbles, is consistent with the expected missing quadratic term in Eq. (|23[) . The coefficient of the linear term for 
spheres, on the other hand, is a factor of 1.6 larger than that for cylinders and is again consistent with the expected 
theoretical factor 2, which is gratifying. Also, our estimates for 5 certainly are compatible with the result of van 
Giessen and Blokhuis (H S = -0.10 ± 0.02. 



7 



IV. RESULTS FROM DENSITY FUNCTIONAL THEORY 



For the one-component Lcnnard- Jones fluid, specified by the interaction potential in Eq. (|17[) . we have performed 
density functional calculations for the metastable bubble and droplets in spherical geometry. Here, we will treat 
the attractive part of the potential in a mean-field fashion which is known to lack quantitative agreement with 
simulations for the phase diagram or the liquid-vapor surface tension. However, the density functional approach 
allows us to study also the limit of large droplet radii to check the validity of the asymptotic expressions ([3]) and 
(|4]), thus complementing our simulation results. 

As usual, the free energy functional of the fluid is split into an ideal part and an excess part, 

F[p] = P d [p]+F^[p] (25) 
with the exact form of the ideal part given by 

F d [p] = Jd 3 rf d (v) = J d 3 rp(r) (ln[p(r)A 3 ] - l) . (26) 

Here, A is the de-Broglie wavelength. The excess part is split into a reference hard-sphere part and an attractive 
part: 

T cx [p] = T hs [p}+ T att [p] . (27) 

For the hard-sphere part we fix the reference hard-sphere diameter to a, for simplicity, and employ fundamental 
measure functionals which are known to be very precise in various circumstances (for recent reviews see Refs. (69l . 

/3JF hs = J dr4>({n[p(r)]}), (28) 

$({n p(r }) = -n ln(l-n 3 + ; + <P ("a) 0A n ^~ • 

Here, $ is a free energy density which is a function of a set of weighted densities {n(r)} = {no, n\, n 2 , n^, ni, n 2 } with 
four scalar and two vector densities. These are related to the density profile p(r) by n a (r) = J dr'p(r') w a (r — r'). 
The weight functions, {w(r)| = {w°, w 1 , w 2 , w 3 , w 1 , w 2 }, depend on the hard sphere radius R = a/2 as follows: 



w 3 = e{R - |r|) , w 2 =S(R- |r|) , w 1 = — , u; 



w 



2 



AttR ' 4ttR 2 ' 



w 2 = ^-M), w-^. (29) 

The reference hard-sphere free energy functional in Eq. ([28)1 is completed upon specification of the function f{n^). 
With the choice <p = 1 we obtain the original Rosenfeld (RF) functional [71(, consistent with the Percus-Yevick 
equation of state. Upon setting 

-t -2n 3 + 3nj-2(l-n3) 2 ln(l-n 3 ) 
3^ (30) 

we obtain the White Bear (WB) functional [Z^lzl], consistent with the quasi-exact Carnahan-Starling equation of 
state. For the attractive part of the free energy functional we employ a mean-field (RPA) approximation motivated 
by the WC A potential separation 74\ : 

F att [p\ = lJ dr J dr' P {v)p{v') w{\v - r'|) , (31) 

W ^ = t > w) • ( 32 ^ 

Here, r m j n = 2 1//6 cr denotes the minimum location of the truncated Lennard- Jones potential u defined in Eq. ()17[) . 



r 2 
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The phase diagram of the model can be calculated easily by evaluating the free energy functional for constant bulk 
densities and employing the Maxwell construction. The critical temperature, T™ f /e ps 1.15 is about 15% too large, 
a defect well-known in mean-field models. The discrepancies in the phase diagram between our mean-field model 
and simulation can be significantly reduced if the reference hard-sphere functional is used to generate a closure 
for the Ornstein-Zernike equation (751 |76| and the reference hard-sphere diameter is appropriately determined 
through a condition on minimal bulk free energy. However, this approach generates a lot more numerical work and 
the results regarding the surface tension of bubbles and droplets remain qualitatively unchanged (see below). 

Planar and spherical surface tensions are obtained by extremizing the grand potential 

Sl[p] = F[p] -nj p(r)dr (33) 
in the appropriate geometry, which leads to the following equation for the equilibrium density profile /o oq (r) 



Po Sp{r) 

where the asymptotic bulk density is given by po an d p = m(Po) is the corresponding chemical potential (p ex 
denotes its excess part). In planar geometry, the grand potential cxtremum is indeed a minimum if the boundary 
conditions for p^ aax (z) are such that the coexisting vapor density p v is enforced asz4 — oo and the corresponding 
coexisting liquid density pi is enforced as z — > oo. The planar surface tension thus becomes 



Ivl 



n'[ P f*™{z)\ - n'w = n'[ P ^\z)] + p ( Pv )l , (35) 



where Q' = SI I A is the aerial grand potential density, and L is the length of our numerical box in ^-direction. 
In spherical geometry, the asymptotic bulk density po = ]hn T -^ 00 has to be chosen either as an oversaturated 
liquid (for a vapor bubble) or an oversaturated vapor (for a liquid bubble). The corresponding extremal point 
of the grand potential is a saddle point which makes the extremization a bit more cumbersome. Simple iteration 
schemes will always diverge ("evaporating" the droplet) after some initial period of convergence. Therefore we have 
chosen to adopt a Picard iteration scheme along with suitable radial shifts of the whole density profile each time 
the Picard iteration starts to diverge. This procedure is repeated until the necessary radial shifts are below the 
resolution of our radial grid which is 0.001 a. From the equilibrium profile Poq h ( r ) the equimolar droplet radius R 
is determined by the condition of no net adsorption. For the associated surface tension it is necessary to introduce 
the asymptotic "inner" density of the drop p\ ^ po which is given by the condition p(p±) = p(po)- (For small 
drops, the actual density in the center is not necessarily equal to p\.) The excess pressure within the drop is defined 
via our mean-field equation of state as Ap = p{p\) — p{po)- Furthermore it is convenient to introduce the excess 
grand potential of the drop by 

4-7T 

AO ^n[pH h (r)} + Y L 3 sph p (36) 

where £ sp h is the radius of our spherical numerical box. This excess grand potential constitutes the nucleation 
barrier at pressure po- The (equimolar) surface tension then becomes 

lvl[R) = 4^{ An + AP ) ' (37) 

The radius R p of the "surface of tension" is defined by the Laplace condition Ap = 2j v i(R p )/ R p . Little algebra 
leads to the relations: 

( 3AS1 \ * 

= \2nAp ) ■ < 38 ' 

- (^y . 

We report numerical results using the WB functional as the reference hard-sphere functional (see Eqs. ([28| and 
(|30| ). Coexistence densities and planar surface tensions for the two investigated temperatures T/T c = 0.68 and 
0.78 are given in Tab. HI Since we are quite far from the critical point, the match in the coexistence densities 



T/T c 


p v a 3 




p„(j 3 p t a 3 7„!cr 2 /e 






(DFT) 


(Simulation) 


0.68 


0.005 


0.79 0.56 


0.010 0.77 0.47 


0.78 


0.010 


0.71 0.39 


0.027 0.71 0.29 



TABLE I: Comparison between DFT and simulation results for the coexistence densities and liquid-vapor surface tension 
at the two investigated temperatures. 

between DFT and simulation is quite reasonable. There is, however, a mismatch in the planar surface tension of 
about 25%. Next we evaluated the surface tensions of bubbles and droplets for radii between 2 and approximately 
200 a. This allows us to test the asymptotic expressions for the Tolman length given in Eqs. © and The 
results for 5(R) = R- R P {R) and ^ v i/lvl{R) - 1 arc shown in Figs. Q7] (for T/T c = 0.78) and[TH] (for T/T c = 0.68). 
One can see that through a linear fit to 8 as a function of 1/R one obtains the Tolman length to a good precision: 
6 = -0.127 ± 0.002 for T/T c = 0.78 and 8 = -0.122 ± 0.002 for T/T c = 0.68. The linear terms in the quadratic 
fits to "f v i/^ v i(R) — 1 (whose modulus should equal 28) are in agreement with these results although here the 
discrepancy between the bubble and droplet results are a bit larger: 8 = —0.129 ± 0.002 for T/T c = 0.78 and 
8 = —0.125 ± 0.008 for T/T c = 0.68. Overall these results for the Tolman length are consistent with the simulation 
results of the previous section, as is the qualitative behavior of "f v i / Jvi(R) — 1- The most notable discrepancy arises 
in the fit coefficient of the quadratic term which is about a factor 4 smaller in DFT, and thus the variation of 
7„;(i?) with R is much less pronounced in DFT than it is in the simulation results. To understand this discrepancy 
partly, note that in the determination of -f v i (R) in DFT (Eq. (|37|) we use the mean-field bulk equation of state in 
the metastable domain in order to subtract the bulk grand potential of the droplet. This can be expected to lead 
to different results compared to the simulations where the free energy density /l {p, T) of a finite system (with box 
length L and density p, exhibiting a stable droplet) is determined and used to subtract the bulk grand potential 
of the droplet. However, one can also expect that the correlations inside the droplet are not adequately captured 
by the mean-field treatment. 

Finally we compare our DFT approach to existing DFT approaches in the literature. All DFT results have 
been obtained usin g a mean-field ansatz for the effect of attractions and differ mainly by the treatment of the 
hard repulsions p5 T l28l - [30l l33l HH, [36|, H(| • Differences in the latter manifest themselves in the absolute numbers 
for the surface tension whereas the results for the reduced quantity ~ivi/lvi{R) are basically unaffected (see e.g. 
Refs. [lHHil). (Differences arise for smaller R upon using different subtraction schemes for the bulk grand potential 
of the droplet). All DFT results for Lennard- Jones type liquids (with differing cut-offs, applied also at various 
temperatures) have so far predicted a negative Tolman length. Here we have confirmed this finding by considering 
also significantly larger droplets and by checking the consistency of the Tolman length extraction from both vapor 
bubbles and liquid droplets. 

V. CONCLUSIONS 

In the present work a computer simulation approach to explore the surface free energy of curved interfaces is 
presented and applied to perform a comparative study of droplets of the minority phase in an unmixed symmetrical 
binary Lennard- Jones mixture and of both droplets and bubbles in the two-phase coexistence region of a simple 
Lennard Jones fluid. In the latter case, both droplets and bubbles of spherical and cylindrical shapes have been 
analyzed with the basic finding of a systematic difference between the curvature-dependent surface free energy of 
droplets and bubbles. Finally, also a density functional treatment of both droplets and bubbles in simple fluids is 
presented that supports the latter observation. 

For the symmetric binary Lennard- Jones mixture, such a difference cannot occur, of course, since an interchange 
of A and B in the two-phase coexistence turns the minority phase into the majority phase, and a droplet becomes a 
bubble, or vice versa. As a consequence, Eq. with a nonzero Tolmann length 8, cannot apply, furthermore the 
dependence of the surface tension 7ab(-R) on the droplet radius of curvature cannot contain any odd term in l/R, 
since R changes sign when a droplet is turned into a bubble. We find that for temperatures far below criticality, 
the simple formula 'Jab(R) = 7ab(oo)/[1 + 2(£ S /R) 2 ] {Eq. (fT6|) } is a very good representation of our data, for radii 
R down to about 3 Lennard-Jones diameters, and the characteristic length £ s is a bit longer than a Lennard- Jones 
diameter (Fig. 
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The main result for the droplets and bubbles in the simple Lennard- Jones fluid (which lacks the trivial particle- 
hole symmetry of the simplistic lattice gas model) is the finding that there is a significant difference between 
7(i?) of droplets and bubbles (Fig. I13|14|15|16[) . This shows that the conventional "capillarity approximation" of 
the standard classical nucleation theory is not strictly valid, since the approximation does not allow for any such 
difference. Also theories such as those of McGraw and Laaksonen [f58j], which imply that the Tolmann length S 
is strictly zero and the difference 7(00) — 7(-R), to leading order, is 1/i? 2 , are hence inconsistent with our results. 
We do find that this difference should be described both by a linear term {S/R) and a quadratic term (£ S /R) 2 
[Eq. fHI ], while the length £ s again is of the order of the Lennard- Jones diameter a, as for the symmetric binary 
Lennard- Jones mixture, the Tolman length 5 being an order of magnitude smaller, 5 ~ —0.16. This order of 
magnitude is compatible with the conclusion of some of the previous work [40L l42j . Qualitatively, these conclusions 
are fully confirmed by the density functional treatment (Sec. IV): note that because of the mean-field character of 
the latter, it does not precisely reproduce the equation of state of the model studied in the simulation and so it 
would be premature to ask for a quantitative "fit" of the simulation results (Fig. [T5l[T6|) by the theory fFigs. [T71[T8|) . 
However, the qualitative similarity is striking, and the order-of-magnitude agreement between the corresponding 
numbers is, of course, gratifying. 

An important consequence of our results also is that deviations from classical nucleation theory for bubble 
nucleation of vapor should be distinctly larger than for droplet nucleation of liquid, at comparable radius R. 
Classical nucleation theory overestimates the nucleation barrier in the range of practical interest. 

Of course, it would be interesting to study the temperature dependence of the lengths 5 and £ s , £ c , and to extend 
our study to other models of fluids. This will be left to future work. 

Acknowledgements: This work was supported in part by the Deutsche Forschungsgcmcinschaft (DFG) through 
the Collaborative Research Centre SFB-TR6 under grants No. TR6/A5 and TR6/N01 and through the Priority 
Program SPP 1296 under grants No. Bi 314/19 and Schi 853/2. S.K.D. is grateful to the Institut fur Physik 
(Mainz) for the hospitality during his extended visits. We are also grateful to the Zentrum fur Datenverarbeitung 
(ZDV) Mainz and the Jiilich Supercomoputer Centre (JSC) for computer time. One of us (K.B.) is grateful to C. 
Dellago, D. Frenkel, G. Jackson and E.A. Miiller for stimulating discussions. 



11 



hi I I I I I I I 



1.4 r 
1.3 r 
1.2 - 



: / 
- / 
1.1 -® 
-I 

1^ 



P 











i i i i i i-i 



6 



^ -. 

Q : 

b- 
\ 

\ - 

\- 



-I 



Q Q(Bjl " I I I I I I I I I I I I I I I I I I I I ij) 







0.2 0.4 0.6 0.8 



FIG. 1: Phase diagram of the symmetric binary (A,B) Lennard- Jones mixture, cf. Eqs. (f7))- (|10[) , in the plane of variables 
T and relative concentration of A particles (xa = Na/N; N = Na + Nb) for fixed N = 6400. The cross shows the critical 
point, as obtained previously [47j |. The horizontal broken line means that phase coexistence is studied at T — 1.0. 
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FIG. 2: Effective free energy /l(xa,T) of finite-size boxes of linear dimension L with L = 10 plotted vs. a; a at T = 1.0, for 
the model of Fig. [1] The estimation of the size-dependent interfacial tension yAB(L) is indicated. 




1/L 

FIG. 3: Extrapolation of 'jab(L) as function of 1/L, cf. Eq. (|12|l . in order to estimate jab(oo). 




FIG. 4: Effective free energy /l(xa,T) plotted vs. natT = 1.0 for L = 10, 12, 16, 26, as indicated. 
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(a) 




FIG. 5: a) Snapshot of a spherical droplet configuration formed by A particles in the background of B-particles (not shown), 
for T — 1.0, L — 24, and xa = 0.15. b) Same as (a), but for a cylindrical droplet, choosing xa = 0.27. Note that our method 
does not at all suppress statistical fluctuations in the size and shape of these droplets, which therefore have spherical or 
cylindrical symmetry on the average only. 
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FIG. 6: Plot of A/j,(Xa) /IcbT vs. Xa at T — 1.0 for L = 10, 12, 16, 26. Data refer to a single run at each size, to illustrate 
the typical noise level (200 Monte Carlo steps per particles have been used for each window of the successive umbrella 
sampling). For the final analysis, 5 such runs were averaged over. 




FIG. 7: Schematic explanation of how the estimation of the functions A/j,l(xa,T) and /l(xa,T) together allows the 
estimation of the concentration difference Aa: and free energy difference A/ due to a droplet. 
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FIG. 8: Plot of Fs/ksT of spherical A-rich droplet at T — 1 for the binary symmetric LJ mixture of Fig.[T] The description in 
terms of the capillarity approximation of classical nucleation theory (CNT) is shown as a broken curve, using 7ab(oo) = 0.722 
as obtained in Fig. [3] The full curve is a superposition of independent simulation results for L — 12, 16,18,20,22 and 24, 
where a running averaging was done using the combined data set. 
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FIG. 9: Plot of 7ab(oo)/7ab(-R) - 1 versus l/R 2 . Here jab(oo) is taken from Fig. El while jab(R) = F s (R)/4:TvR 2 is 
estimated using Eq. (|14|) . Ideally the estimates obtained from different values of L should superimpose on a single curve. 
The scatter between the curves for different values of L is due to residual statistical errors. The thin straight line is a fit 
function giving 7as(oo)/7.4b(-R) — 1 ~ 2.2/ R 2 . 
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FIG. 10: Effective free energy density /i,(p)/fcsT of the single-component Lennard- Jones fluid at T = 0.78T C plotted vs. 
density p for 3 values of L, as indicated in the figure. 
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FIG. 11: Plot of pl{p) /ksT for the one-component Lennard- Jones fluid as a function of p, at T = 0.78T C , for 3 values of L, 
as indicated. 
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FIG. 13: Plots of Fs/ksT of spherical droplets and bubbles for the one-component LJ fluid, at T = 0.78T C , as a function 
of sphere radius R. The capillarity approximation (CNT), Fs/ksT = 4irR 2 y v i(oo) is included, using the estimate of 7(00) 
from Fig. [12] 




FIG. 14: Same as Fig. 1131 but for cylindrical droplets and bubbles. Note that, here the j/-axis corresponds to the surface 
free energy per unit height of the cylinder. 
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FIG. 15: Plots of y v i(R)/y v l(oo) - 1 vs. 1/R for spherical droplets and bubbles for the LJ fluid at (a) T = 0.78T C and (b) 
T = 0.68T C . Fits to the functional forms (|20|l are included. 
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FIG. 16: Same as Fig. 1151 but for cylindrical droplets and bubbles. 
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FIG. 17: (a) Radius-dependent Tolman length S(R) = R — R p for bubbles and droplets with a corresponding linear fit in 
the range 1/R £ (0,0.1). (b) Surface tension ratio y v l/jvl(R) — 1 as function of the equimolar radius R for bubbles and 
droplets with a corresponding quadratic fit in the range 1/R G (0, 0.1). All results are for the temperature T = 0.78T C . 




FIG. 18: (a) Radius-dependent Tolman length S(R) = R — R p for bubbles and droplets with a corresponding linear fit in 
the range 1/R £ (0,0.1). (b) Surface tension ratio jt,l/jvl(R) — 1 as function of the equimolar radius R for bubbles and 
droplets with a corresponding quadratic fit in the range 1/R G (0,0.1). All results are for the temperature T = 0.68T C . 
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